

#include <functional>
#include <iostream>

typedef std::function<double(double)> mathfunc;

double trapezoid( mathfunc f, double a, double b )
{
  return (b - a) * ((f(a) + f(b)) / 2);
}

double simpson( mathfunc f, double a, double b )
{
  return ((b - a) / 6) * (f(a) + 4 * f((a+b)/2) + f(b));
}

double integral( mathfunc f, double a, double b, double delta = 0.001 )
{
  if( b - a < delta )
    return (b - a) * ((f(a) + f(b)) / 2);

  auto m((a + b) / 2);
  return integral(f, a, m, delta) + integral(f, m, b, delta);
}

typedef std::function<double(mathfunc,double,double)> method;

double integral( method r, mathfunc f, double a, double b, double delta = 0.001 )
{
  if( b - a < delta )
    return r(f, a, b);

  auto m((a + b) / 2);
  return integral(r, f, a, m, delta) + integral(r, f, m, b, delta);
}


double q( double x )
{
  return x * x;
}

/*
 * 
 */
int main( int argc, char** argv )
{
  std::cout << integral(q, -1, 1) << std::endl;
  std::cout << integral(simpson, q, -1, 1) << std::endl;
  return 0;
}

